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Numerical study of boundary layer interaction 
with shocks - method and code validation 

By N. A. Adams 
1. Motivation and objectives 

A major problem in modeling of turbulent supersonic flows is the correct assess- 
ment of viscous-inviscid interaction problems. Of particular interest is the interac- 
tion of boundary layers with shocks. Present turbulence models give in most cases 
unsatisfactory results in the region of rapid distortion and in the separation region 
(if one is present) in particular with regard to mean flow profiles and turbulence 
quantities (cf. Kline et a/., 1981). 

Recent direct numerical simulations (DNS) at moderate supersonic Mach num- 
bers of boundary layers without interaction show that the effect of compressibility 
in those cases is rather small (Guo k Adams, 1994). Even at those Mach num- 
bers, however, compressibility can have a significant effect in case viscous-inviscid 
interaction is present. 

Compression corner flows are of great practical interest since they appear to be 
omnipresent in aeronautical configurations (aircraft fuselage, fuselage-wing junc- 
tion, engine inlet, etc.). On the other hand, they also give rise to a particularly 
interesting combination of phenomena, which all are more or less confined to a 
relatively narrow region about the corner. First, the turbulence in the oncom- 
ing boundary layer responds to a rapid distortion. This is a generalization of the 
problem of isotropic turbulence interacting with normal shocks (e.g. Lee, 1993) to 
anisotropic turbulence with inhomogeneous mean shear. Second, for large enough 
Mach numbers and deflection angles there is a shock-induced unsteady separation. 
The separation bubble is contained by a detached curved shear layer, and fluc- 
tuations in this shear layer, subject to high strain, are strongly amplified. Some 
experiments report evidence for Gortler-like vortices in the detached shear layer 
(Smits k Muck, 1987). Third, there is unsteady shock motion, which is suspected 
to be triggered by bursting events in the oncoming turbulent boundary layer (An- 
dreopoulos k Muck, 1987). And finally, the shock is generated at the wall. This 
is favorable for direct numerical simulations since it relieves the need to accurately 
introduce a shock at the outer boundaries. 

The objective of the present work is the direct numerical simulation of shock 
boundary layer interaction. This report summarizes the first phase during which 
a numerical method suitable for this problem has been developed and a computer 
code has been written and tested. 

2* Accomplishments 

The first part of this work focuses on the development of a new type of spatial 
discretization scheme which combines the spectral-like wave- represent at ion of sym- 
metric compact finite difference schemes (Lele, 1992) with a suppression of aliasing 
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errors by implicit dissipation at non-resolved wavenumbers and increased robustness 
against generation of spurious waves at the boundaries. In an additional step the 
scheme is made shock capturing by switching locally to an essentially non-oscillatory 
(ENO) scheme near discontinuities (in case of weak solutions). In this report the 
numerical method is merely outlined; for details the reader is referred to Adams & 
Shariff (1994). 

The newly developed spatial scheme is used to discretize the reduced hyperbolic 
part of the Navier-Stokes equations while for the parabolic part a symmetric com- 
pact finite difference scheme is used (cf. Adams, 1993). This report summarizes 
the status after implementation and validation of the according computer code. 

The performance of the present working code lies in between about 200 MFLOPS 
and about 140/js per grid point and time step (pure 3rd order ENO, 2D) and about 
500 MFLOPS and roughly 20 fis per grid point and time step (pure compact-FD, 
3D) on a single CRAY Y-MP C90 CPU. These relations reflect the increased number 
of logical instructions required for the non-linear ENO scheme. 

2.1 Numerical method 

Generalizing the formulation of compact finite- difference schemes (Lele, 1992), 
a family of centered upwind-biased compact schemes of 5th order is introduced. 
Numerical dissipation is used to suppress unresolved wavenumbers while an accurate 
representation of the dispersion relation for resolved wavenumbers is required. This 
requirement is formulated as a constrained optimization problem. A parameter 
study allows for the generation of a whole family of locally optimal schemes, one 
member of which, called P455/1, we are using in this study. The general formula 
of the schemes is 

5Z a ^j+M “ a vfj+v • (1) 

V=-Vl 

In this equation fj = f(xj) is the grid function, the derivative of order a of which 
is searched for. For scheme P455/1 it is a = 1, /i/ = /i r = v\ = v r = 2 for interior 
schemes, /i/ = v\ = 1, fi r = 2, v r = 3 for the left next -to- boundary scheme and 
m = vi = 0, fx r = 2, v r = 4 for the left boundary scheme (accordingly for the 
schemes at the right boundary). The coefficients and a u are determined from 
order conditions and the abovementioned optimization problem. For the numerical 
values see Adams & Shariff (1994). 

The discrete derivative operator for a semi- discretized scalar advection equation 
on a strip consists of interior and boundary schemes whose frequency responses are 
shown in Fig. 1. It should be noted that the dispersion, Fig. la, is well approximated 
up to wavenumbers larger than 2, significantly increasing the resolved- wavenumber 
domain when compared with non-compact finite-difference schemes (e.g. that of 
Rai & Moin, 1993). The numerical dissipation is shown in Fig. lb. 

Moreover, the compact formulation allows for stable high-order boundary closures 
avoiding the order-drop at the boundary as in Rai & Moin (1993). Note that it 
has been shown by Gustafsson (1975) that for hyperbolic equations the boundary 
closures must not be of order less than (r — 1) to maintain a global order r of 
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FIGURE 1. Frequency response u r of scheme P455/1, (a) real part, (b) imaginary 
part. scheme at left boundary; scheme next to left boundary; in- 
terior scheme; scheme at right boundary; scheme next to right boundary. 

£ is the wavenumber normalized with the grid spacing. 


the semidiscretization. The linear stability of this scheme is analyzed in Adams k 
Shariff (1994). 

Two different approaches of making an underlying compact scheme shock captur- 
ing (nonlinearly stable) have been pursued initially. The first followed the approach 
of Cockburn k Shu (1994), which required some minor modifications to be appli- 
cable for our general type of upwind compact schemes. Numerical tests, however, 
showed an unsatisfactory shock resolution. A Gibbs-like phenomenon could not be 
suppressed satisfactorily without a significant smearing of the shock. The second 
follows the guideline of Hou k Le Floch (1994), who gave evidence that a noncon- 
servative scheme converges to a weak solution (so a solution exists) if in the neigh- 
borhood of discontinuities a conservative scheme is used. In our case scheme P455/1 
is used in the smooth regions, while a 5th order ENO scheme in finite-difference 
form (Shu k Osher, 1989) using Roe fluxes on local characteristics with entropy-fix 
by switching to a local Lax- Friedrichs flux formulation is used near discontinuities. 
Without being able to give theoretical evidence, numerical tests reported in Adams 
k Shariff (1994) demonstrate that the hybrid scheme possesses the ENO property. 

For integration in time two different 3rd order Runge-Kutta methods are used. A 
pure ENO scheme as spatial discretization is combined with a TVD Runge Kutta 
scheme (Shu, 1988). The hybrid scheme is combined with a low-storage explicit 
Runge Kutta scheme (Williamson, 1980, case 7). The admissible time step is calcu- 
lated from an estimated bound for the discrete convective and the discrete viscous 
operator using a CFL (Courant-Friedrichs-Lewy) condition. 

2.2 Discretization of the compressible conservation equations 

The fundamental equations are the volume- specific conservation equations for 
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mass momentum and energy 


dv _ dF OG m 

dt dx dy dz + 


(i) 


The flux vectors F, G, and H are given, for example, in Anderson et ai( 1984). 
The viscosity /i is calculated from Sutherland’s law. Z is a distributed forcing 
which can be used to cancel the residual of a given basic flow (this is done in the 
calculations of sections 2.4.1 and 2.4.2, otherwise Z — 0). For a discussion of this 
term see, for instance, Adams (1993). 

The fundamental systems of PDE are discretized in a method of lines manner, i.e. 
the spatial derivatives of the fluxes are approximated by a spatial discretization. The 
system thus becomes an ODE in t in terms of the grid point values and is projected 
forward in time with a Runge-Kutta integration method. 

The viscous fluxes are discretized in conservative form by a symmetric P346- 
scheme used in Adams (1993). Note that although a conservative discretization of 
a linear heat equation with this scheme would be asymptotically unstable (Adams, 
1993), this kind of discretization remains bounded in the present case with an 
upwind discretization of the convective terms. Since a conservative discretization 
has significantly less operations, it has been given the preference over the non- 
conservative discretization used in Adams (1993). 

The discretization of the convective terms has been a major objective. For the 
details the reader is referred to Adams & Shariff (1994). Here, it is summarized in 
a few words. At each t the fluxes calculated from the instantaneous solution are 
processed by a discontinuity detector algorithm, and cells are marked for treatment 
by the ENO scheme. In the present code this is done in each index space component, 
and a whole plane is marked if at least one cell of this plane contains a transition. 
This is favorable for vectorization but has the disadvantage that the ENO scheme 
may be used in smooth regions, too. For the marked cells the fluxes are projected on 
the local characteristics by a transformation with the left-eigenvector modal matrix 
of the cell’s Roe matrix. The flux derivatives at the cell faces (i.e. the nodes) 
are reconstructed from the cell-centered numerical fluxes on local characteristics 
obtained with a 5th order ENO scheme using a Roe-flux formulation with a local 
Lax-Friedrichs flux as entropy fix. The numerical fluxes (after application of the 
ENO procedure) are then projected back onto the computational space basis by a 
transformation with the corresponding right-eigenvector modal matrix. 

The flux derivatives in the smooth regions are calculated with the positive biased 
and the negative biased compact schemes. This is done by projecting the fluxes 
with the respective right-hand side matrices of the schemes onto a local average, 
where the entries at the grid points at which the solution is to be taken from the 
ENO procedure are replaced with the already known flux derivatives. The left- 
hand side is modified accordingly by setting the respective submatrix to unity. In 
the reconstruction step the so obtained linear equation systems (pentadiagonal) are 
solved for the positive and negative biased fluxes. The upwinding is done by choos- 
ing the flux derivative in upwind direction according to the direction of the local 
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characteristics at each grid point. In case of a vanishing eigenvalue, the arithmetic 
mean is taken as is for the flux derivatives obtained by ENO since they are already 
upwinded. The advantage of this non-building block (though straight forward) flux 
splitting is that the flux-approximation is smooth up to the order of the scheme at 
sonic points. This would not be the case in general if split fluxes would be used. 

2.S Outflow boundary treatment 

A major concern of direct numerical simulation methods for spatially evolving 
convective problems is the correct formulation at the outflow boundary. Since this 
boundary is artificial a prescription of the correct conditions would require the 
knowledge about the solution at this plane. Several approaches are now in more or 
less standard use. For incompressible flows there is the relaminarization method of 
Kloker et al.{ 1993). In compressible flow Pruett et a/.(1994) use a buffer domain 
approach. Guo et a/.(1994) suggested a much simpler method with the same effi- 
ciency as the buffer domain approach. In this work we adopt the latter, though in 
a different formulation, to account for the non- perturbation form of the equations 
solved here since the basic flow residual is not compensated by a forcing term as 
in Pruett et a/., 1994, and Guo et a/., 1994. Simpler approaches, for instance non- 
reflecting conditions (e.g. Poinsot & Lele, 1992), treat viscous wave forms within 
the boundary layer improperly and give rise to spurious reflected waves. 

The basic idea follows closely the concept of a sponge layer according to Israeli 

Orszag (1981). From a simple model equation as 


dv dv 
dt dx 


a(x)v 


( 2 ) 


it is seen that cr(x) in the last term of the right-hand side has the character of a 
Newtonian cooling coefficent. 

We choose the following damping function (cf. Israeli &; Orszag, 1981) 


<t(x)=A 3 (N s + 1)(N 9 +2) 


(x - x 3 ) n ’(L x - a:) 
( Lx ~ x s ) N *+ 2 


(3) 


for x 3 < x < L x . It has the properties: (i) j^ z a(x)dx = A s and (ii) cr(L x ) = 0. 
The implementation is done readily by adding a term 


Z 3 = —cr(x)(U — Uo) 

for x 9 < x < L x to the fundamental equations (1). Herein U denotes the vector of 
conservative variables and Uq = f7|* = o. From property (ii) of a it is clear that the 
fundamental equations (1) are recovered at x = L x where a perfect non-reflecting 
condition is prescribed (Thompson, 1987). 

2.4 Code validation 

The computer code is validated using two kind of tests. The first is typical for 
DNS codes where a known mean flow is enforced and the correct representation of a 
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TABLE 1. Flow parameters. 


Case 

Moo 

•Res, 

T./Too 

a 

P 

L z 

A 

2.5 

12773 

2.041 

0.26 

0.43 

15 

B 

4.5 

10000 

4.38 

2.25 

0 



linear eigensolution is required. The second is typical for CFD codes where a mean 
flow is to be found. 

In the first case we check for the efficiency and correct implementation of deriva- 
tive and integration routines by marching eigensolutions of a temporally evolving 
boundary layer (streamwise periodic) in time. The global growth rate (defined from 
the mode energy Eq. (4)) is required to match the real part of the eigenvalue. Also, 
since the amplitude function is only a function of the wall normal coordinate z , fre- 
quency and growth rate at a given z should be close to the imaginary and real part 
of the eigenvalue locally in z. Inflow and outflow boundary conditions are examined 
by considering spatial instability in a parallel boundary layer. The spatial growth 
rate is required to be sufficiently accurately represented and the outflow boundary 
condition should have a limited upstream effect. 

In the second test the steady problem of a shock impinging on a flat plate bound- 
ary layer is investigated. Though none of the methods used here is suitable to obtain 
steady state solutions efficiently, evidence is given that the (quasi time dependent) 
solution marches toward a steady state. After a reasonable number of iterations 
the computations are halted and compared with reference data. 

2.^.1 Eigensolutions of temporally evolving boundary layer 

As test cases we choose cases B and C from Adams (1993, section 6.1). For 
completeness the flow parameters are given in table 1. A combined al gebraic- si n/i 
mapping is used in z; see Adams, 1993, section 6.1 for parameters. Boundary 
conditions are periodic in x and y , non- reflecting at z = L Z} and isothermal no-slip 
at the wall with the wall temperature set to the adiabatic wall temperature (Mack, 
1984). 

The unperturbed mean flow is calculated from the compressible similarity equa- 
tions (Stewartson, 1964) with a shooting method (Adams, 1993). An unstable 
eigenmode with streamwise and spanwise wavenumbers a and /?, respectively, as 
given in table 1 calculated from a spectral solution method (Simen, 1993) (used 
earlier in Adams, 1993) is superimposed. For case A the eigensolution is an oblique 
vortical mode, for case B a two dimensional Mack mode (mixed vortical/acoustic). 
In both cases the initial amplitude is A — 1CT 4 . The grid spacing in x and y is 
uniform, the box dimensions are chosen consistently with the wavenumbers such 
that the eigenmodes become (±1,±1) and (±1,0) Fourier modes, respectively. 

For the calculations scheme P455/1 is used. Since the solution is smooth the 
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(a) (b) 


FIGURE 2. Global growth rate, (a) test case A, (10 x 10 x 50), (20 x 

20 x 50), (20 x 20 x 100), (40 x 40 x 100), — - linear theory, (b) 

test case B, (10 x 50), (20 x 50), (20 x 100), (40 x 100), 

— — ’ linear theory. 

ENO scheme is inactive. We define a modal energy as 

I* L x 

E(k z , k y \t) = / p(z)uj (&x i hy ! t)uj (k x , ky , t)dz (4) 

Jo 

Uj stands for the three velocity components (sum over j). 

Fig. 2 shows the growth-rate obtained from the modal energy for different dis- 
cretizations for 1100 time steps each. The initial transient is an effect of different 
solution methods and different meshes used for the initial eigenmode and for the 
DNS. 

Figs. 3 and 4 show the local growth rate and frequency across the boundary layer 
for test cases A and B, respectively, making use of the fact that the eigenfunction 
maintain their shape in a parallel boundary layer. The improvement from a dis- 
cretization 10 x 50 to 20 x 100 is mainly due to the refinement in s which allows 
for a better resolution of the region around the critical layer which is also a region 
of high curvature of the mean-flow profiles. 

2-4-2 Eigensolutions of spatially evolving boundary layer 

Since the j/-discretization remains unchanged compared to the test cases in sec- 
tion 2.4.1, we can restrict ourselves here to 2D-problems. The main concern is 
the correct formulation of the inflow boundary condition and the efficiency of the 
outflow- boundary condition given in section 2.3. The test case is the spatially evolv- 
ing equivalent of case B of section 2.4.1. Also, since it has been found in section 
2.4.1 that N z = 51 is about sufficient to resolve the eigenfunctions in 2 , we focus 
on the effect of changing N x and the outflow boundary condition parameters. At 
the inflow all primitive variables axe prescribed (as function of t) according to the 
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(a) (b) 

FIGURE 3. Local growth rate and frequency test case A. (a) growth rate u (b) 

frequency uv. />, (20 x 50), u, (20 x 20 x 50), p, (40 x 40 x 100), 

u, (40 x 40 x 100), 1 1 1 linear theory. 




(*) (b) 

FIGURE 4. Local growth rate and frequency, test case B, (a) growth rate u;,. 

(b) frequency u > T . p, (20), u, (20 x 50), p, (40 x 100), u, 

(40 x 100), linear theory. 


well-posedness requirement (Oliger & Sundstrom, 1978); at the outflow the sponge 
layer of section 2.3 together with a non-reflecting condition at x = L x is used. As 
validation-test we check for the correct spatial growth rates in terms of primitive 
variables at an arbitrary position z. The growth rates are obtained from the coeffi- 
cients of the solution’s Fourier transform in y and t. These coefficients assume the 
form 

Oyi{x,z) = A{z)e ,az (5) 

where A{z) is the complex amplitude function of the linear eigensolution. From two 
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(a) ( b ) 

FIGURE 5. Local growth rate, N x = 81, A s = 4, N s = 3, z = 0.48, St = 4tt/uj. (a) 
sponge layer and non-reflecting boundary condition, (b) non-reflecting boundary 
condition only. p, u, — linear theory. 


successive ^-stations the complex a can then be obtained by 


a = -i 2 , ( 6 ) 

x 2 1 x i 

so that the growth rate is — Im(a) and the frequency Re(a). 

Fig. 5 shows the effect of the sponge layer. From the experience with a symmetric 
compact scheme (Guo et a/., 1994) we choose the sponge layer thickness to be one 
wavelength of the eigenmode. The minimum necessary sponge-layer thickness is 
problem dependent and is therefore not further investigated here. Spurious reflected 
waves penetrate into the computational domain further upstream in the case of a 
pure non- reflecting boundary condition, Fig. 5b, than in the case of a non-reflecting 
boundary condition plus a sponge layer, Fig. 5a. Though this effect is not so 
pronounced for the low amplitude linear perturbation in the present case, it is 
expected to be stronger for large amplitude turbulent fluctuations. Note that non- 
reflecting boundary conditions are not consistent with the evolution of a boundary 
layer eigensolution. This reflects in large oscillations near the outflow boundary of 
a sensitive measure as the local growth rate, Fig. 5b. 

From Fig. 6 it is evident that increasing the order N s of the cooling function 
polynomial in Eq. (3), (making the damping more localized) or increasing the cool- 
ing intensity A s has no noticeable effect on the growth rate or the extent of the 
valid domain. Thus it is expected that an increased damping is able to extinguish 
stronger turbulent fluctuation efficiently without increasing the invalid part of the 
domain. 

Fig. 7a shows the same case as Fig. 5a, but the integration time is increased 
to seven periods so that the wave can travel about two times the width of the 
computational box. Obviously the region of upstream influence of the boundary 
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FIGURE 6. Local growth rate, N x = 81, St = 4n/u;. (a) A s = 4, N s = 7. (b) 
A s = 10, A T S = 3. p, a, 1 linear theory. 



(a) (b) 

FIGURE 7. Local growth rate, (a) N x = 81, N z = 51, 6i = 147 rfu. (b) N x = 161, 
iV z — 101, = 4n/uj. /?, u, 1 linear theory. 


condition did not increase, so that the shorter integration times used before lead 
to correct conclusions about the upstream extent of the invalid part of the domain. 
In Fig. 7b both streamwise and vertical number of grid points are doubled, which 
results in a more accurate approximation of the growth rate as expected. 

Two remarks are in order. Firstly we note that the setup used in the above 
test calculations is quite severe, allowing for only 4 wavelengths A of the primary 
wave in z, which finally results in about 2.5A for the extent of the valid part of 
the domain. Second, inflow transients are apparently unavoidable (as the initial 
transients were in the preceding section) since the inflow perturbation is taken 
from a linear stability solution obtained with a different method and on a different 
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TABLE 2. Flow parameters. 


quantity 

value 

comment 


221. 6K 

estimated 

Moo 

2 


P*o 0 

9178. 79Pa 

estimated 

Pr 

0.72 


K 

1.4 


R 

287.03 


Plo 

1.449 • 10 -5 kg/m/s 


s* 

110.4K 


Reto 

140000 


Re (l 

450000 


Re 6l 

909.9053 


6* 

1.5309 • 10 -4 m 

at £ch from sim. sol. 

K 

4.7325 • 10 -4 m 

99.9% thickness 

io 

82.7559 


L z 

415 


L z 

200 



mesh. Other DNS results (e.g. Pruett et al, 1994, Guo et al., 1994) confirm this 
observation. 

It can be concluded that the outflow boundary treatment has no upstream effect 
further than roughly 1.5A upstream of L x for ( L x — == A. The findings in this 

section resemble the experiences made by Guo et a/.(1994) with a different spatial 
discretization and a different cooling term. Pruett et a/.(1994) report that their 
buffer domain approach spoils an upstream region of about 2 A. The sponge layer 
approach thus allows for a much simpler formulation with a comparable performance 
of a well-tuned buffer-domain approach. 

2>4 -3 Laminar boundary layer interacting with an impinging shock 

The following test example is a standard case for the validation of steady state 
Navier-Stokes solvers. Experimental data are provided by Hakkinen et a/.(1959). 
Although the experimental evaluation is limited and the comparison suffers from 
incompletely reported flow parameters, the experiment’s favorable feature is that 
the flow is laminar, although it should be noted that the test section extends into 
the region where the laminar flat plate boundary layer is unstable. 

An extensive numerical investigation of this particular problem has been done 
by Katzer (1989). We emphasize here that for the results presented in this section 
time-accurate and low-dissipation methods have been used. The computations have 
thus been halted before a true steady state has been reached. The flow parameters 
are given in table 2 (reference length is 6*, dimensional quantities are marked with 
a star). 
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82 100 


498 


FIGURE 8. Quasi-Schlieren plot (intensity proportional to norm of density gradi- 
ent). 


Fig. 8 shows a quasi- Schlieren plot (merely the norm of the density gradient) 
when the computations were halted. A shock, introduced at the inflow boundary, 
impinges at (3 = 32.6° on the laminar boundary layer along a flat plate. At the 
boundary layer edge it is reflected as a decompression wave. The shock-induced 
boundary layer separation gives rise to a compression wave in front of the separation 
region and a compression wave behind. 

As initial condition we take the solution of the inviscid shock-reflection problem 
outside of the boundary layer, while near the wall a boundary layer from a similarity 
solution is given. The shock is to impinge on the plate at x sh = 325.2041 for the 
inviscid problem. 

As boundary conditions we fix at the inflow the initial condition for all primitive 
variables (giving the correct number of 5 conditions for the Navier-Stokes equations 
according to Oliger & Sundstrom, 1978). At the outflow we prescribe perfectly 
non-reflecting boundary conditions (Thompson, 1987), and no viscous conditions 
are imposed. In fact the imposition of weak viscous boundary conditions in terms 
of derivatives of the viscous fluxes (Poinsot & Lele, 1992) was found to have no effect, 
while the imposition of boundary conditions in terms of stresses (Dutt, 1988) re- 
sulted in an outflow boundary-layer due to the inconsistency between boundary con- 
dition and solution. At the wall we prescribe a no-slip adiabatic condition. At the 
upper boundary all flow variables are prescribed corresponding to the state behind 
the impinging shock. Non-reflecting conditions at the upper boundary (Thompson, 
1987) were found to give rise to a viscous (heat-equation like) instability emerging 
from the corner between inflow and upper boundary. The reason for that behavior 
is that the presence of the shock close to the upper boundary at the inflow leads 
to non-negligible viscous terms near the edge of the computational box, and inflow 
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Figure 9. Residual. Method A: 3rd order ENO; method B: 5th order hybrid. 
p, pu, pw , — - E. 

and boundary conditions which are well posed only in the inviscid limit become 
invalid. 

We perform a calculation with N x = 151 and N y = 100. First, the meaningless 
transient caused by the initial condition is spanned by marching a 3rd order ENO 
scheme 5500 iterations in time (method A). Note that since we use a time accurate 
solution procedure and since the ENO-stencil changes in time, the residual does not 
reach machine zero but merely the truncation error of the spatial discretization. 
This has been confirmed by continuing the calculation with method A for 5500 
additional iterations (not shown). After the residual settled down the computation 
is continued for another 5500 iterations with the hybrid scheme (method B). Since 
the Reynolds number is small, the shock is resolvable by the scheme and the ENO- 
scheme is only active in z in the outer 4.95% of the domain (in the average over 
all iterations), where the grid spacing is wide (shock detector parameter settings 
: fi x = 1 , f} z = 0.05). Fig. 9 shows the evolution of the residual for the four 
conservative variables. 

Figs. 10 and 11 show pressure contours after 5500 iterations with method A and 
method B, respectively. The hybrid scheme appears to represent the reflected com- 
pression and decompression waves more accurately than the pure ENO scheme. This 
suggests that the internal phenomena of the boundary layer are better represented. 




352 


N. A. Adams 



FIGURE 10. Pressure contours after 5500 iterations, method A ( min = 0.16, 
max = 0.26, inc = 0.002). 



FIGURE 11. Pressure contours after 5500 additional iterations, method B (min = 
0.16, max = 0.26, inc = 0.002). 

This finding is confirmed if we compare numerical and experimental results. 
Fig. 12 shows surface pressure and skin friction. In both cases method B gives 
a better representation. Note that the experimental values downstream of the sepa- 
ration region appear to be affected by the particular method of measuring the skin 
friction with a pressure probe (Hakkinen, 1959). The accurate numerical solutions 
of Katzer (1989) show that the experimental values in this region are somewhat too 
large. 
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(a) * 

FIGURE 12. Surface pressure (a) 


after A after B. 



x 


(V 

skin friction (b). symbols experiment, 



Figure 13. Velocity profiles, (a) after A (b) after B. The number-pairs on 
the curves give Re X ezp (first entry) and f?e Ieip /fie Inum (second entry), symbols 
experiment, computation. 


The correspondence between experimentally and numerically obtained velocity 
profiles is reasonable for method B and less satisfactory for method A, figure 13. 

Two remarks are in order. First, as can be seen from Fig. 9, the residual is 
still decreasing when the computation with method B was halted, though with a 
shallow slope as expected for a time-accurate method. Thus the results should not 
be considered as converged. Second, the residual level after the computation with 
method A was halted does not decrease if the computation is continued with method 
A. The further decrease of the residual and the improvement of the solution is due 
to method B. 

We finally note that for the present 2D calculations the code performance of the 
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pure 3rd order ENO scheme was 228 MFLOPS and 143 ns per grid point and time 
step, while the 5th order hybrid scheme performed at 354 MFLOPS and 23 n$ per 
grid point and time step (CRAY fpp and eft ^ 7 optimization only, single processor 
Y-MP C90). 

3. Future plans 

In the next step the computational code will be extended to generalized coor- 
dinates (in (x,r)). This requires the characteristic transformation routines to be 
adapted and an extension of the flux calculation routines. 

The final objective is to simulate a compression corner flow according to a suitable 
experiment. Experiments presently considered are those of Smits & Muck (1987) 
(having the advantage of relatively detailed turbulence data, the disadvantage of a 
quite large Reynolds number Re 6i = 78000) of Ardonceau et a/, (lacking the details 
of the previous one but having a smaller Reynolds number of roughly Re$ 2 = 7700), 
and of Zheltovodov et a/. (at about Re Sl = 9720 for the incoming boundary layer 
with the problem that mean flow and turbulence data have been obtained in different 
wind tunnels). The Mach number range of those experiments is between 2.25 and 
3, the turning angles go from 8° to 25°. 
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